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Glossary and Notation 

Discrete variational mechanics: A formulation of mechanics in discrete-time that is based on a 
discrete analogue of Hamilton's principle, which states that the system takes a trajectory for 
. which the action integral is stationary. 

' Geometric integrator: A numerical method for obtaining numerical solutions of differential equa- 

. tions that preserves geometric properties of the continuous flow, such as symplecticity, momentum 

CO ' preservation, and the structure of the configuration space. 

yf^ . Lie group: A differentiable manifold with a group structure where the composition is differentiable. 

' The corresponding Lie algebra is the tangent space to the Lie group based at the identity element. 

. Symplectic: A map is said to be symplectic if given any initial volume in phase space, the sum of the 

' signed projected volumes onto each position-momentum subspace is invariant under the map. One 

consequence of symplecticity is that the map is volume-preserving as well. 

• 1— I 

^ ' 1. Definition of the Subject and Its' Importance 

H : . 

5t 1 Discrete control systems, as considered here, refer to the control theory of discrete-time Lagrangian or 

Hamiltonian systems. These discrete-time models are based on a discrete variational principle, and are part 
of the broader field of geometric integration. Geometric integrators are numerical integration methods that 
preserve geometric properties of continuous systems, such as conservation of the symplectic form, momentum, 
and energy. They also guarantee that the discrete flow remains on the manifold on which the continuous 
system evolves, an important property in the case of rigid-body dynamics. 

In nonlinear control, one typically relies on differential geometric and dynamical systems techniques to 
prove properties such as stability, controllability, and optimality. More generally, the geometric structure of 
such systems plays a critical role in the nonlinear analysis of the corresponding control problems. Despite 
the critical role of geometry and mechanics in the analysis of nonlinear control systems, nonlinear control 
algorithms have typically been implemented using numerical schemes that ignore the underlying geometry. 

The field of discrete control system aims to address this deficiency by restricting the approximation to 
choice of a discrete-time model, and developing an associated control theory that does not introduce any 
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additional approximation. In particular, this involves the construction of a control theory for discrete-time 
models based on geometric integrators that yields numerical implementations of nonlinear and geometric 
control algorithms that preserve the crucial underlying geometric structure. 



2. Introduction 



The dynamics of Lagrangian and Hamiltonian systems have unique geometric properties; the Hamiltonian 
flow is symplectic, the total energy is conserved in the absence of non-conservative forces, and the momentum 
map associated with a symmetry of the system is preserved. Many interesting dynamics evolve on a non- 
Euclidean space. For example, the configuration space of a spherical pendulum is the two-sphere, and 
the configuration space of a rigid body attitude dynamics has a Lie group structure, namely the special 
orthogonal group. These geometric features determine the qualitative behavior of the system, and serve as 
a basis for theoretical study. 

Geometric numerical integrators are numerical integration algorithms that preserve stru ctures of the 



contin uous dynamics such as invariants, symplecticity, and the configuration manifold (see iHairer et al 



( 20061 )). The exact geometric properties of the discrete flow not only generate improved qualitative behavior, 
but also provide accurate and efficient numerical techniques. In this article, we view a geometric integrator 
as an intrinsically discrete dynamical system, instead of concentrating on the numerical approximation of a 
continuous trajectory. 

Numerical integration method s that preserve the simplecticity of a Hamiltonian system have been studied 
(see Leimkuhler and Reich ( 2004 ): Sanz-Serna ( 19921 )'). Coefficients of a Runge-Kutta method are carefully to 
chosen to satisfy a simplecticity criterion and order conditions to obtain a symplectic Runge-Kutta method. 
However, it can be difficult to construct such integrators, and it is not guaranteed that other invariants of 
the system, such as a momentum map, are preserved. Alternatively, variational integrators are constructed 
by discretizing Hamilton 's principle, rather than di scretizing the continuous Euler-Lagrange equation (see 
Marsden and Wesd (|200l[ ): lMoser and Veselovl (|l99lh ). The resulting integrators have the desirable property 
that they are symplectic and momentum preserving, and they exhibit good energy behavior for exponentially 
long times. Lie group methods are nume rical integrators that preserve the Lie group structure of the 
configuration space (see llserles et al. ( 2000( )). Recently, these two approaches have been unified to obtain 
Lie group variational integrators that preserve the geometric properties of the dynamics as well as the Lie 
group structure o f the configu r ation space without th e use of local charts, reprojection, or constraints (see 
Lee et all (|2007al ) : Ilbo^ (|2004l ) : iMarsden et all (|l999f )). 

Optimal control problems involve finding a control input such that a certain optimality objective is 
achieved under prescribed constraints. An optimal control problem that minimizes a performance index is 
described by a set of differential equations, which can be derived using Pontryagin's maximum principle. 
Discrete optimal control problems involve finding a control input for a discrete dynamic system such that 
an optimality objective is achieved with prescribed constraints. Optimality conditions are derived from 
the discrete equations of motion, described by a set of discrete equations. This approach is in contrast 
to traditional techniques where a discretization appears at the last stage to solve the optimality condition 
numerically. Discrete mechanics and optimal control approac hes determin e optim al control inputs and 
trajectories more accurately with less computational load (see Junge et al.l (l2005l)). Combin e d with an 
i ndirec t optim ization technique, they are substantially more efficient (see iHussein et ah ( 2006 ): Lee et al 
(l2005al.l200d)V 



The geometric approach to mechanics can provide as the theoretical basis of innovative control method- 
ologies in geometric control theory. For example, these techniques allow the attitude of satellites to be 
controlled using changes in its shape, as opposed to chemical propulsion. While the geometric structure of 
mechanical systems plays a critical role in the construction of geometric control algorithms, these algorithms 
have typically been implemented using numerical schemes that ignore the underlying geometry. By applying 
geometric control algorithms to discrete mechanics that preserve geometric properties, we obtain exact nu- 
merical implementation of the geometric control theory. In particular, the method of controlled Lagrangian 
systems is based on the idea of adopting a feedback control to realize a modification of either the potential 
energy or the kinetic energy, referred to as potential shaping or kinetic shaping, respectively. These ideas 
are applied to const ruct a real-t i me d igital feedback controller that stabilizes the inverted equilibrium of the 
cart-pendulum (see Bloch et al. (|2005. . .2006) ). 
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Figure 1 . Procedures to derive continuous and discrete equations of motion 



In this article, we will survey discrete Lagrangian and Hamiltonian mechanics, and their applications to 
optimal control and feedback control theory. 

3. Discrete Lagrangian and Hamiltonian Mechanics 

Mechanics studies the dynamics of physical bodies acting under forces and potential fields. In Lagrangian 
mechanics, the trajectory of the object is derived by finding the path that minimizes the integral of a 
Lagrangian over time, called the action integral. In many classical problems, the Lagrangian is chosen as the 
difference between kinetic energy and potential energy. The Legendre transformation provides an alternative 
description of mechanical systems, referred to as Hamiltonian mechanics. 

Discrete Lagrangian and Hamiltonian mechanics has been developed by reformulating the theorems and 
the procedures of Lagran gian and Hamiltonian mechanics in a discrete time setting (see, for example, 
Marsden and Wesd (1200 ih l. Therefore, discrete mechanics has a parallel structure with the mechanics de- 



scribed in continuous time, as summarized in Figure[T]for Lagrangian mechanics. In this section, we describe 
discrete Lagrangian mechanics in more detail, and we derive discrete Euler-Lagrange equations for several 
mechanical systems. 

Consider a mechanical system on a configuration space Q, which is the space of possible positions. The 
Lagrangian depends on the position and velocity, which are elements of the tangent bundle to Q, denoted 
TQ. Let L : TQ — > R be the Lagrangian of the system. The discrete Lagrangian, : Q x Q ^ R is an 
approximation to the exact discrete Lagrangian, 

(1) ir'*(go,9l)- r L{qoi{t),qoi{t))dt, 

Jo 

where qoi(O) — Qo^ Qaiih) = <Zi, and qoi{t) satisfies the Euler-Lagrange equation in the time interval (0, h). 
A discrete action sum ®d ■ Q^^^ R, analogous to the action integral, is given by 

(2) ©d(go, gi, ■ • ■ , 9Ar) = X! Ld{qk,qk+i)- 

The discrete Hamilton's principle states that 

6(3d = 
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for any 5qk, which yields the discrete Eulei — Lagrange (DEL) equation, 
(3) D2La{qk-i,qk) + DiLd{qk,qk+i) = 0. 

This yields a discrete Lagrangian flow map {qk-i,qk) ilkiqk+i)- The discrete Legendre transformation, 
which from a pair of positions {qo, qi) gives a position-momentum pair {qo,Po) — (qo, —DiLd{qo, qi) provides 
a discrete Hamiltonian flow map in terms of momenta. 

The discrete equations of motion, referred to as variational integrators, inherit the geometric properties 
of the continuous system; they are symplectic, and they preserve any momentum maps associated with 
symmetries as the discrete Noether's theorem is satisfied. They exhibit good total energy behavior for 
exponentially long time periods. 

Many interesting Lagrangian and Hamiltonian systems, such as rigid bodies evolve on a Lie group. Lie 
group variational integrators preserve the nonlinear structure of the Lie g roup confi gurations as well as 



hgurationg 

geometric properties of the continuous dynamics (see iMarsden et al.l (|l999t ) and iLeoJ (j2004h ). The basic 



idea for all Lie group methods is to express the update map for the group elements in terms of the group 
operation, 

(4) 51 = 50 /o, 

where go, gi G are configuration variables in a Lie group G, and /o G G is the discrete update represented 
by a right group operation on 50. Since the group element is updated by a group operation, the group 
structure is preserved automatically without need of parameterizations, constraints, or re-projection. In the 
Lie group variational integrator, the expression for the fiow map is obtained from the discrete variational 
principle on a Lie group, the same procedure presented in Figure[T] But, the infinitesimal variation of a Lie 
group element must be carefully expressed to respect the structure of the Lie group. For example, it can be 
expressed in terms of the exponential map as 

d 

g exp £77 = 577, 

=0 

for a Lie algebra element ij G Q. This approach has been appli ed to the rotation group SO (3) and t o the 



tor a Lie algebra element ?7 t 0. inis approacn Has been appli ea to tne rotatiori group bU(dl and t o tne 
special Euclidean group SE(3) for dynamics of rigid bodies (see Lee et al. ( 2005b l. lLee et al" ( 2007a[ ). and 



iLee et al ] (|2007bh 'l. Generalizations to arbitrary Lie groups gives the generalized discrete Euler-Poincare 



(DEP) equation, 

(5) T;L/„ • DMgo, fo) - Ad}„ • (r,*L/, • DMg,,h)) + TtL^, ■ D^L,{g^,h) = 0, 

for a discrete Lagrangian on a Lie group, : G x G ^ M. Here Lf : G G denotes the left translation 
map given hy Lfg = fg for f,g€G, TgLf : TgG TfgG is the tangential map for the left translatio n, and 
Adg : ^ is the adjoint map. A dual map is denoted by a superscript * (see Marsden and Ratiul ([l999) 



for detailed definitions). 

We illustrate the properties of discrete mechanics using several mechanical systems, namely a mass-spring 
system, a planar pendulum, a spherical pendulum, and a rigid body. 

Example 1 (Mass-spring System). 

Consider a mass-spring system, defined by a rigid body that moves along a straight frictionless slot, and is 
attached to a linear spring. 

Continuous equation of motion: The configuration space is Q = M, and the Lagrangian L : ffi. x M ^ M is 
given by 

(6) L{q,q) ^^mf -^Kq'^, 

where g G M is the displacement of the body measured from the point where the spring exerts no force. 
The mass of the body and the spring constant are denoted by m, k € M, respectively. The Euler-Lagrange 
equation yields the continuous equation of motion. 

(7) mq + Kq = 0. 

Discrete equation of motion: Let /i > be a discrete time step, and a subscript k denotes the A:-th discrete 
variable at t = kh. The discrete Lagrangian : R x M — > K is an approximation of the integral of the 
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(a) Mass-spring system 
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(b) Planar pendulum 



Figure 2. Computed total energy (RK45: blue, dotted, VI: red, solid) 
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along the solution of ([7|) over a time step. Here, we choose the following discrete 



(8) 



Ldiqk,qk+i) ^ hL 



qk + qk+i qk+i - qk 



qk) 



hn 



qk+if 



Direct application of the discrete Euler-Lagrange equation to this discrete Lagrangian yields the discrete 
equations of motion. We develop the discrete equation of motion using the discrete Hamilton's principle 
to illustrate the principles more explicitly. Let : ^ M be the discrete action sum defined as 

C5rf = X^fclo^ Lii{qk,qk+i), which approximates the action integral. The infinitesimal variation of the action 
sum can be written as 

AT-l 



60d = ^ Sqk+i \ ^{qk+i - qk) ~^{qk + qk+i) \ + 5qk 

k=0 



Since Sqo = Sqj\j 



h '"'^"^ 4 
0, the summation index can be rewritten as 



'-{qk+i - qk) 



hn 



(qk + qk+i) 



N-l f 



-((Jfc+i - 2% + qk-i) ^(ft+i + 29fc + qk-i) 



From discrete Hamilton's principle, 5<S>a 



(9) 



Thus, the discrete equation of motion is given by 

'-{qk+i + 2gA; + qk-i) = 0. 



for any Sqk- 

m, hn 

-rilk+i - ^qk + qk-i) + -r 
h 4 

For a given (<7fe-i, qk), we solve the above equation to obtain qk+i- This yields a discrete fiow map (gfe-i, qk) ^ 
{qk, qk+i), and this process is repeated. The discrete Legendre transformation provides the discrete equation 
of motion in terms of the velocity as 

4m 



1 + ] qk+i = /ifjfe + ( 1 - ) qk 



hn hn 
-Ik - —Qk+i- 



(10) 

2m- 2m ^ 

For a given {qk, qk), we compute qk+i and (jk+i by (|10p and pip , respectively. This yields a discrete flow map 
{qk,qk) ^ {qk+i,qk+i)- it can be shown that this variational integrator has second-order accuracy, which 
follows from the fact that the discrete action sum is a second-order approximation of the action integral. 

Numerical example: We compare computational properties of the discrete equations of motion given by 
([TO)l and pTjl with a 4(5)-th order variable step size Runge-Kutta method. We choose m = 1 kg, k ~ 1 kg/s^ 
so that the natural frequency is 1 rad/s. The initial conditions are qo — \/2m, 9o = 0; and the total energy is 
E = 1 Nm. The simulation time is 2007r sec, and the step-size h = 0.035 of the discrete equations of motion 



is chosen such that the CPU times are the same for both methods. Figure 2(a) shows the computed total 
energy. The variational integrator preserves the total energy well. The mean variation is 2.7327 x 10""'^^ Nm. 
But, there is a notable dissipation of the computed total energy for the Runge-Kutta method 



6 



TAEYOUNG LEE, MELVIN LEOK, AND N. HARRIS MCCLAMROCH 



Example 2 (Planar Pendulum). 

A planar pendulum is a mass particle connected to a frictionless, one degree-of-freedom pivot by a rigid 
massless link under a uniform gravitational potential. The configuration space is the one-sphere S"'^ = 
jgeM^I ||g||=l}. While it is common to parameterize the one-sphere by an angle, we develop parameter- 
free equations of motion in the special orthogonal group S0(2), which is a group of 2 x 2 orthogonal matrices 
with determinant 1, i.e. S0(2) = {i? € M2x2 | ^ j^^^^ det[i?] = l}. S0(2) is diffeomorphic to the one- 
sphere. It is also possible to develop global equations of motion on the one-sphere directly, as shown in 
the next example, but here we focus on the special orthogonal group in order to illustrate the key steps to 
develop a Lie group variational integrator. 

We first exploit the basic structures of the Lie group SO (2). Define a hat map which maps a scalar 
to a 2 X 2 skew-symmetric matrix f2 as 

"0 -n 
n 



n = 



The 2x2 skew-symmetric matrices are the Lie algebra so (2). Using the hat map, we identify so (2) with R. 
An inner product on so (2) can be induced from the inner product on K as 

any 17i,J72 G The matrix exponential is a local diffeomorphism from so(2) to S0(2) given by 

cos ft — sin ri 
sin n cos il 



exp Cl 



The kinematics equation for R e SO (2) can be written in terms of a Lie algebra element as 



(12) 



R = Rn. 



Continuous equations of motion: The Lagrangian for a planar pendulum L 
written as 



S0(2) X so(2) 



can be 



(13) 



L(i?, Cl) = Knfn'^ + mgle. 



Re2 — —mP 



mgle2 Re2, 



where the constant 5 € R is the gravitational acceleration. The mass and the length of the pendulum are 
denoted by m, / S R, respectively. The second expression is used to define a discrete Lagrangian later. We 
choose the bases of the inertial frame and the body-fixed frame such that the unit vector along the gravity 
direction in the inertial frame, and the unit vector along the pendulum axis in the body-fixed frame are 
represented by the same vector 62 = [0; 1] £ K^. Thus, for example, the hanging attitude is represented by 
R = I2x2- Here, the rotation matrix R € SO (2) represents the linear transformation from a representation 
of a vector in the body-fixed frame to the inertial frame. 

Since the special orthogonal group is not a linear vector space, the expression for the variation should be 
carefully chosen. The infinitesimal variation of a rotation matrix R G SO (2) can be written in terms of its 
Lie algebra element as 



(14) 



SR 



d 



Rcxpef] ~ Rrj, 



where rj d R so that ^ G so (2). The infinitesimal variation of the angular velocity is induced from this 
expression and (fT2l) as 

5R^R + R^5R 



(15) 



5VL 



-t)Q. + VLfj + f] = r). 



Define the action integral to be © = L{R, Cl) dt. The infinitesimal variation of the action integral is 
obtained by using (|14p and (jlSp . Hamilton's principle yields the following continuous equations of motion. 

(16) 

(17) R^Rn. 



n + je^Rei = 0, 



If we parameterize the rotation matrix as R 
(18) 



cos f — sm f 
sin 9 cos 9 



these equations are equivalent to 



= 0. 
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Discrete equations of motion: We develop a Lie group variational integrator on SO (2). Similar to 
define Fk £ S0(2) such that 

(19) ^fc+i — RkFk- 

Thus, Fk — R^Rk+i represents the relative update between two integration steps. If we find Fk G S0(2), the 
orthogonal structure is preserved through (|19p since multiplication of orthogonal matrices is also orthogonal. 
This is a key idea of Lie group variational integrators. 

Define the discrete Lagrangian Ld : SO (2) x SO (2) ^ R to be 

Ld{Rk,Fk) = -^"^^^{^k - l2x2,Fk - 12x2) + ^mgle2Rke2 + ^mgle2 Rk+ie2, 

(20) = -^mlhv[l2x2 - Fk] + hmgle'^Rke2 + ^mgle2 Rk+ie2, 

which is obtained by an approximation hClk — R]^{Rk+i — Rk) — Fk ~ I2x2, applied to the continuous 
Lagrangian given by ((T3|). 

As for the continuous time case, expressions for the infinitesimal variations should be carefully chosen. 
The infinitesimal variation of a rotation matrix is the same as (jl4p . namely 

(21) SRk = RkVk, 

for r]k € K, and the constrained variation of Fk is obtained from (fT9|) as 

(22) 6Fk = SRlRk+i + Rl5Rk+i = -i%Fk + Fkf]k+i = Fkim+i - F^fjkFk) = Fk{f]k+i - Vk), 

where we use the fact that FfjF^ — fj for any F E S0(2) and fj E so(2). 

Define an action sum i3d : 80(2)^+^ ^ M as (S^ = J2k=o Ld{Rk,Fk). Using ^ and (HH), the variation 
of the action sum is written as 

^-1/1 I 
^^d^Y.\ ^"^^'(^''-1 ~ ^^-1) - ^{Fk - Fl) - hmgle^RkCi, r% 

k=l ^ 

From discrete Hamilton's principle, (525^ = for any f/k- Thus, we obtain the Lie group variational integrator 
on S0(2) as 

(23) (Fk ~ Fj) - {Fk+i - Fj+i) - ^e^R^ie, = 0, 

(24) Rk+i = RkFk. 

For a given {Rk,Fk) and Rk+i = RkFk, (|23| is solved to find Fk+i. This yields a discrete map {Rk,Fk) 1-^ 
{Rk+i, Fk+i). If we parameterize the rotation matrices R and F with 9 and A0 and if we assume that 
A6 <^ 1, these equations are equivalent to 

^{0k+i - 29k + Ok) + Y sin^fe = 0. 
The discrete version of the Legendre transformation provides the discrete Hamiltonian map as follows. 

(25) Fk - = 2hn - l^efRTei, 

(26) -Rfc+i — RkFk, 

(27) f)fc+i = ilfc - ■^e^i?feei - ■^e^i?fc+iei. 

For a given {Rk,^k), we solve ((25|) to obtain Fk- Using this, {Rk+i,^k+i) is obtained from ((26l) and (|27p . 
This yields a discrete map (Rk, ^k) ^ {Rk+i, ^^fc+i)- 

Numerical example: We compare the computational properties of the discrete equations of motion given 
by (P5)) - ([77|) with a 4(5)-th order variable step size Runge-Kutta method. We choose m = 1 kg, I = 9.81 m. 
The initial conditions are = 7r/2rad, = 0, and the total energy is i? = ONm. The simulation time is 
1000 sec, and the step-size h — 0.03 of the discrete equations of motion is chosen such that the CPU times 
are identical. Figure 2(b) shows the computed total energy for both methods. The variational integrator 
preserves the total energy well. There is no drift in the computed total energy, and the mean variation is 
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1.0835 X 10 ^ Nm. But, there is a notable dissipation of the computed total energy for the Runge-Kutta 
method. Note that the computed total energy would further decrease as the simulation time increases. 

Example 3 (Spherical Pendulum). 

A spherical pendulum is a mass particle connected to a frictionless, two degree-of-freedom pivot by a rigid 
massless link. The mass particle acts under a uniform gravitational potential. The configuration space is the 
two-sphere = {9 G | \\q\\ = l}- It is common to parameterize the two-sphere by two angles, but this 
description of the spherical pendulum has a singularity. Any trajectory near to singularity causes numerical 
ill-conditioning. Furthermore, this leads to complicated expressions involving trigonometric functions. 

Here we develop equations of motion for a spherical pendulum using the global structure of the two-sphere 
without parameterization. In the previous example, we develop equations of motion for a planar pendulum 
using the fact that the one-sphere is diffeomorphic to the special orthogonal group S0(2). But, the 
two-sphere is not diffeomorphic to a Lie group. Instead, there exists a natural Lie group action on the two- 
sphere. That is the 3-dimensional special orthogonal group SO (3), a group of 3 x 3 orthogonal matrices with 
determinant 1, i.e. S0(3) ^ {R £ M^xs | ^ /g^g, det[i?] = l}. The special orthogonal group S0(3) acts 
on the two-sphere in a transitive way; for any qi,q2 G there exists a i? € S0(3) such that q2 = Rqi- 
Thus, the discrete update for the two-sphere can be expressed in terms of a rotation matrix as This is 

a key idea to develop a discrete equations of motion for a spherical pendulum. 

Continuous equations of motion: Let g € be a unit vector from the pivot point to the point mass. The 
Lagrangian for a spherical pendulum can be written as 

(28) L{q,q) = ^ml'^q- q + mgles- q, 

where the gravity direction is assumed to be 63 = [0; 0; 1] € K'^. The mass and the length of the pendulum 
are denoted by m, Z € M, respectively. The infinitesimal variation of the unit vector q can be written in terms 
of the vector cross product as 

(29) Sq = ^xq, 

where ^ e M"^ is constrained to be orthogonal to the unit vector, i.e. ^ ■ q — 0. Using this expression for the 
infinitesimal variation, Hamilton's principle yields the following continuous equations of motion. 

(30) q+{q-q)q+^{qx{qxe3))^0. 

Since q — u x q for some angular velocity G M'^ satisfying uj ■ q = 0, this can also be equivalently written as 

(31) u;-|gxe3 = 0, 

(32) q ^ UJ X q. 

These are global equations of motion for a spherical pendulum; these are much simpler than the equations 
expressed in term of angles, and they have no singularity. 

Discrete equations of motion: We develop a variational integrator for the spherical pendulum defined on 
§^ . Since the special orthogonal group acts on the two-sphere transitively, we can define the discrete update 
map for the unit vector as 

(33) qk+i = Fkqk 

for Fk e SO (3). Define a discrete Lagrangian : §^ x ^ K to be 

Ld{qk,qk+i) = -^mf{qk+i - qk) ■ (qk+i - qk) + ^mgles ■ qk + ^mgles ■ qk+i- 
The variation of qk is the same as (1^^ . namely 



(34) 



5qk ^ £,k 'x qk 
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^1 ^ "^0 50 100 150 200 50 100 150 200 

(a) Trajectory of pendulum (b) Total energy (c) Unit length error ■ 5fc — 1| 



Figure 3. Numerical simulation of a spherical pendulum (RK45: blue, dotted, VI: red, solid) 



for e M'^ with a constraint • = 0. Using this discrete Lagrangian and the expression for the variation, 
discrete Hamilton's principle yields the following discrete equations of motion for a spherical pendulum. 



(35) 
(36) 



+ X 63 I X gffe + II 

hg 



1/2 



hg 



UJk+l = UJk + -^Qk X 63 + -^Qk+l X 63, 



Since an explicit solution for Fk € SO (3) can be obtained in this case, the rotation matrix Fk does not 
appear in the equations of motion. This variational integrator on exactly preserves the unit length of qk , 
the constraint qk ■ oJk = 0, and the third component of the angular velocity uJk ■ 63 which is conserved since 
gravity exerts no moment along the gravity direction 63. 

Numerical example: We compare the computational properties of the discrete equations of motion given 
by ([55)1 and ([55]) with a 4(5)-th order variable step size Runge-Kutta method for PT|) and ([5^ . We choose 

m = 1kg, / = 9.81m. The initial conditions are go — ["^jOj^]; ^0 = 0.1[\/3, 0, 3] rad/sec, and the total 
energy is E = — 0.44Nm. The simulation time is 200 sec, and the step-size h = 0.05 of the discrete equations 
of motion is chosen such that the CPU times are identical. Figure [3] shows the computed total energy 
and the unit length errors. The variational integrator preserves the total energy and the structure of the 
two-sphere well. The mean total energy deviation is 1.5460 x 10~^Nm, and the mean unit length error 
is 3.2476 X 10"^^. But, there is a notable dissipation of the computed total energy for the Runge-Kutta 
method. The Runge-Kutta method also fails to preserve the structure of the two-sphere. The mean unit 
length error is 1.0164 x 10^^. 



Example 4 (Rigid Body in a Potential Field). 

Consider a rigid body under a potential field that is dependent on the position and the attitude of the 
body. The configuration space is the special Euclidean group, which is a semi-direct product of the special 
orthogonal group and Euclidean space, i.e. SE(3) = S0(3) ©R^^ 

Continuous equations of motion: The eq uations of motion for a rigid body can be developed either from 
Hamilton's principle (see Lee et al. ( 2007a[) ) in a similar way as Example [21 or directly from the generalized 
discrete Euler-Poincare equation given at ([5])- Here, we summarize the results. Let m e M and J e R'^^'^ be 
the mass and the moment of inertia matrix of a rigid body. For {R,x) e SE(3), the linear transformation 
from the body-fixed frame to the inertial frame is denoted by the rotation matrix R e SO (3), and the 
position of the mass center in the inertial frame is denoted by a vector x e R'^. The vectors 51, w g R^^ are the 
angular velocity in the body-fixed frame, and the translational velocity in the inertial frame, respectively. 
Suppose that the rigid body acts under a configuration-dependent potential U : SE(3) M. The continuous 
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400 500 



(a) Trajectory of rigid body (b) Total energy (c) Rotation matrix error ||/ — 

Figure 4. Numerical simulation of a dumbbell rigid body (LGVI: red, solid, RK45 with 
rotation matrices: blue, dash-dotted, RK45 with quaternions: black, dotted) 



equations of motion for the rigid body can be written as 

(37) R = 

(38) i = V, 

(39) JVl + VLx JVL = A4, 

dU 

(40) mil = - — , 

ox 

where the hat map " is an isomorphism from R'^ to 3 x 3 skew-symmetric matrices so (3), defined such that 
xy = X X y ior any G M'^. The moment due to the potential M e R'^ is obtained by the following 
relationship. 

, , ' dU^ rdU 

The matrix ^ G R'^^^ is defined such that — for i,j G {1, 2, 3}, where the i,j-th element of a 

matrix is denoted by 

Discrete equations of motion: The corresponding discrete equations of motion are given by 

(42) hJQk + —Mk = FkJd - JdFj, 

(43) ^fe+i — RkFk, 

dUk 



(44) Xk+i = Xk + hvk - 



2m dxk ' 



(45) Jilk+i - F^jrik + ^F^Mk + ^Mk+u 

(46) mvk+i = mvk - - -7: , 

2 dxk 2 dxk+i 

where Jd G M.^^^ is a non-standard moment of inertia matrix defined as Jd — ^ti^[J]l3x3 — J- For a 
given (i?fc, Xfc, fife, Wfc), we solve the implicit equation ([1^ to find Fk G SO (3). Then, the configuration 
at the next step Rk+i,Xk+i is obtained by and (UH), and the moment and force Mk+i, g^^^^ can 
be computed. Velocities ilk+i,Vk+i are obtained from pS)) and (pS)) . This defines a discrete fiow map, 
(i?fc, Xfc, ffe) 1-^ (i?fe+i, Xfc+i, f2fc+i, Wfe+i), and this process can be repeated. This Lie group variational 
integrator on SE(3) can be ge neralized to multiple rigid bodies acting under their mutual gravitational 
potential (see lLee et al. I (l2007al) '). 

Numerical example: We compare the computational properties of the discrete equations of motion given 
by (|42 p - ((46|) with a 4(5)-th order variable step size Runge-Kutta method for ([37|) -([40 |l . In addition, we 
compute the attitude dynamics using quaternions on the unit three-sphere S"^. The attitude kinematics 
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equation p7p is rewritten in terms of quaternions, and the corresponding equations are integrated by the 
same Runge-Kutta method. 

We choose a dumbbeh spacecraft, that is two spheres connected by a rigid massless rod, acting under a 
central gravitational potential. The resulting system is a restricted full two body problem. The dumbbell 
spacecraft model has an analytic expression for the gravitational potential, resulting in a nontrivial coupling 
between the attitude dynamics and the orbital dynamics. 



As shown in Figure 4(a) the initial conditions are chosen such that the resulting motion is a near-circular 
orbit combined with a rotational motion. Figure |4(b)| and Figure 4(c) show the computed total energy 
and the orthogonality errors of the rotation matrix. The Lie group variational integrator preserves the total 
energy and the Lie group structure of S0(3). The mean total energy deviation is 2.5983 x 10^'', and the mean 
orthogonality error is 1.8553 x 10~^^. But, there is a notable dissipation of the computed total energy and 
the orthogonality error for the Runge-Kutta method. The mean orthogonality errors for the Runge-Kutta 
method are 0.0287 and 0.0753, respectively, using kinematics equation with rotation matrices, and using 
the kinematics equation with quaternions. Thus, the attitude of the rigid body is not accurately computed 
for Runge-Kutta methods. It is interesting to see that the Runge-Kutta method with quaternions, which 
is generally assumed to have better computational properties than the kinematics equation with rotation 
matrices, has larger total energy error and orthogonality error. Since the unit length of the quaternion vector 
is not preserved in the numerical computations, orthogonality errors arise when converted to a rotation 
matrix. This suggests that it is critical to preserve the structure of S0(3) in order to study the global 
characteristics of the rigid body dynamics. 

The dynamics of a rigid body is characterized by a Hamiltonian system on a Lie group. The Lie group 
variational integrator has a desirable property that both the symplectic structure and the Lie group structure 
of the rigid body dynamics are preserved concurrently. More explicitly, the computational properties of the 
Lie group variational integrator is compared with a symplectic integrator that does not pr eserve the Lie grou p 
structure, and a Lie group method that does not preserve the symplectic structure (see Lee et al. ( 2007bl )). 
It is shown that the Lie group variational integrator has a substantial superiority in terms of numerical 
accuracy and efficiency. Due to these computational advantages, the Lie group variational integrator has 
been used to study the dynamics of the binary near-Earth asteroid 66391 (1 999 KW4) in joint wo rk between 
the University of Michigan and the Jet Propulsion Laboratory, NASA f see IScheeres et aL ( 2006f) ). 



4. Optimal Control of Discrete Lagrangian and Hamiltonian System 

Optimal control problems involve finding a control input such that a certain optimality objective is 
achieved under prescribed constraints. An optimal control problem that minimizes a performance index is 
described by a set of differential equations, which can be derived using Pontryagin's maximum principle. 
The equations of motion for a system are constrained by Lagrange multipliers, and necessary conditions for 
optimality is obtained by the calculus of variation. The solution for the corresponding two point boundary 
value problem provides the optimal control input. Alternatively, a sub-optimal control law is obtained by 
approximating the control input history with finite data points. 

Discrete optimal control problems involve finding a control input for a given system described by discrete 
Lagrangian and Hamiltonian mechanics. The control inputs are parameterized by their values at each discrete 
time step, and the discrete equations of motion are derived from the discrete Lagrange-d'Alembert principle 
(|Kane et all (|2000l )'l. 



JV-l N-l 

S X! Ld{qk,qk+l) = X! ^^dilk, Qk+l) ■ Sqk + F^{qk,qk+i) ' 5qk+l\, 
k=0 k=a 

which modifies the discrete Hamilton's principle by taking into account the virtual work of the external 
forces. Discrete optimal control is in contrast to traditional techniques such as collocation, wherein the 
continuous equations of motion are imposed as constraints at a set of collocation points, since this approach 
induces constraints on the configuration at each discrete timestep. 

Any optimal control algorithm can be applied to discrete Lagrangian or Hamiltonian system. For an 
indirect method, our approach to a discrete optimal control problem can be considered as a multiple stage 
variational problem. The discrete equations of motion are derived by the discrete variational principle. 
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The corresponding variational integrators are imposed as dynamic constraints to be satisfied by using La- 
grange multipliers, and necessary conditions for optimality, expressed as discrete equations on multipliers, 
are obtained from a variational principle. For a direct method, control inputs can be optimized by using 
parameter optimization tools such as a sequential quadratic programming. The discrete optimal control can 
be characterized by discretizing the optimal control problem from the problem formulation stage. 

This method has substantial computational advantage to find an optimal control law. As discussed in the 
previous section, the discrete dynamics are more faithful to the continuous equations of motion, and conse- 
quently more accurate solutions to the optimal control problem are obtained. The external control inputs 
break the Lagrangian and Hamiltonian system structure. For example, the total energy is not conserved for 
a controlled mechanical system. But, the computational superiority of the discrete mechanics still holds for 
controlled systems. It has been shown that the discrete dynamics is more reliable even f or controlled system 



as it computes the energy dissipation rate of controlled systems more accurately (see iMarsden and West 



(j2001l )). For example, this feature is extremely important in computing accurate optimal trajectories for 
long term spacecraft attitude maneuvers using low energy control inputs. 

The discrete dynamics does not only provide an accurate optimal control input, but also enables us 
to find it efficiently. For the indirect optimal control approach, optimal solutions are usually sensitive to 
a small variation of multipliers. This causes difficulties, such as numerical ill-conditioning, to solve the 
necessary conditions for optimality expressed as a two point boundary value problem. Sensitivity derivatives 
along the discrete necessary conditions do not have numerical dissipation caused by conventional numerical 
integration schemes. Thus, they are numerically more robust, and the necessary conditions can be solved 
computationally efficiently. For the direct optimal control approach, optimal control inputs can be obtained 
by using larger discrete step size, which requires less computational load. 

We illustrate the basis properties of the discrete optimal control using optimal control problems for the 
spherical pendulum and the rigid body model presented in the previous section. 

Example 5 (Optimal Control of a Spherical Pendulum). 

We study an optimal control problem for the a spherical pendulum described in Example [3l We assume 
that an external control moment u G M"^ acts on the pendulum. Control inputs are parameterized by their 
values at each time step, and the discrete equations of motion are modified to include the effects of the 
external control inputs by using discrete Lagrange-d'Alembert principle. Since the component of the control 
moment that is parallel to the direction along the pendulum has no effect, we parameterize the control input 
as Uk = qk X Wk for Wk € K^. 

The objective is to transfer the pendulum from a given initial configuration (qq^luq) to a desired configu- 
ration {q'^jOj'^) during a fixed maneuver time N, while minimizing the square of the weighted I2 norm of the 
control moments. 

^ h rp J-^ h , 



mmj = ^ -ujuk = X! 2 ^ ^''^^(^fe ^ 

k=0 k=0 



We solve this optimal control problem by using a direct numerical optimization method. The terminal 
boundary condition is imposed as an equality constraint, and the 3A^ control input parameters {u!k}^^Q 
are numerically optimized using sequential quadratic pro gramming. This me thod is referred to as a DMOC 
(Discrete Mechanics and Optimal Control) approach (see Junge et al. ( 20051 ) V 



Figure [5] shows a optimal solution transferring the spherical pendulum from a hanging configuration given 
by (go,^^o) = (63,03x1) to an inverted configuration {q'^,uj'^) = (—63,63x1) during 1 second. The time step 
size is ft = 0.033. Experiment has shown that the DMOC approach can compute optimal solutions using 
larger step size than typical Runge-Kutta methods, and consequently, it requires less computational load. 
In this case, using a second-order accurate Runge-Kutta method, the same optimization code fails while 
giving error messages of inaccurate and singular gradient computations. It is presumed that the unit length 
errors of the Runge-Kutta method, shown in Example [31 cause numerical instabilities for a finite-difference 
gradient computations required for the sequential quadratic programming algorithm. 

Example 6 (Optimal Control of a Rigid Body in a Potential Field). 

We st udy an opti mal control problem of a rigid body using a dumbbell spacecraft model described in Example 
H (see iLee et all (jioOG, ) for detail). We assume that external control forces u-^ g M'^, and control moment 
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Figure 5. Optimal control of a spherical pendulum 



G M.^ act on the dumbbell spacecraft. Control inputs are parameterized by their values at each time step, 
and the Lie group variational integrators are modified to include the effects of the external control inputs by 
using discrete Lagrange-d'Alembert principle. 

The objective is to transfer the dumbbell from a given initial configuration (Rq,xq,Qq,vq) to a desired 
configuration {R"^ , x'^ , ^l'^ , v'^) during a fixed maneuver time N, while minimizing the square of the I2 norm 
of the control inputs. 



min : 



k=0 



where Wf, Wm S ffi are symmetric positive definite matrices. Here we use a modified version of the discrete 
equations of motion with first order accuracy, as it yields a compact form for the necessary conditions. 

Necessary conditions for optimality: We solve this optimal control problem by using an indirect opti- 
mization method, where necessary conditions for optimality are derived using variational arguments, and a 
solution of the corresponding two-point boundary value problem provides the optimal control. This approach 
is common in the optimal control literatures; here the optimal control problem is discretized at the problem 
formulation level using the Lie group variational integrator presented in Section [3l 

N-l , , 



fe=0 



+ AJ.'^ {-Xk+i + Xk + hvk} + Xl'^ { -mvk+i + mvk - h 



2,T 
k 



dUk+1 



+ Xl'^ (logm(Ffc - RlRk+i)) + Xk^ {-J^k 



fc+i 



dxk+ 
Jilk + h (Mk+i 



1)}- 



where A^, A^, A^, € are Lagrange multipliers. The matrix logarithm is denoted by logm : SO (3) — > so (3) 
and the vee map V : so (3) is the inverse of the hat map introduced in Example |31 The logarithm 

form of (|43p is used, and the constraint (|42p is considered implicitly using constrained variations. Using 
similar expressions for the variation of the rotation matrix and the angular velocity given in (jl4[) and (jl5p . 
the infinitesimal variation can be written as 



k=l 



[Ai.;A^;A^;At] e and zk e 

= [\og-sn{Rl5RkY ■,5xk,5nk,5vk]. The matrix G 



where A^ = 
given by Zk 

Thus, necessary conditions for optimality are given by 

4+1 = -wf'xl 

rn T/f/"! \4 

Afc = Afc^iAfe+1 



represents the infinitesimal variation of {Rk,Xk,^k,Vk), 
12x12 jg defined in terms of {Rk, Xk,^k,Vk)- 



(47) 
(48) 
(49) 
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(a) Convergence rate (b) Orbital radius change (c) Inclination angle change 

Figure 6. Optimal control of a rigid body 



together with the discrete equations of motion and the boundary conditions. 

Computational approach: Necessary conditions for optimality are expressed in terms of a two point bound- 
ary problem. This problem is to find the optimal discrete flow, multipliers, and control inputs to satisfy the 
equations of motion, optimality conditions, multiplier equations, and boundary conditions simultaneously. 
We use a neighboring extremal method fsee lBrvson and Hoi ( 19751 )'). A nominal solution satisfying all of the 
necessary conditions except the boundary conditions is chosen. The unspecified initial multiplier is updated 
by successive linearization so as to satisfy the specified terminal boundary conditions in the limit. This is 
also referred to as a shooting method. The main advantage of the neighboring extremal method is that the 
number of iteration variables is small. 

The difficulty is that the extremal solutions are sensitive to small changes in the unspecified initial mul- 
tiplier values. The nonlinearities also make it hard to construct an accurate estimate of sensitivity, thereby 
resulting in numerical ill-conditioning. Therefore, it is important to compute the sensitivities accurately 
to apply the neighboring extremal method. Here the optimality conditions (|47|) and ((48|) are substituted 
into the equations of motion and the multiplier equations, which are linearized to obtain sensitivity deriva- 
tives of an optimal solution with respect to boundary conditions. Using this sensitivity, an initial guess of 
the unspecified initial conditions is iterated to satisfy the specified terminal conditions in the limit. Any 
type of Newton iteration can be app l ied. W e use a line search with backtracking algorithm, referred to as 
Newton-Armijo iteration (see Kelley ( 19951 )). 

Figure [6] shows optimized maneuvers, where a dumbbell spacecraft on a reference circular orbit is trans- 



ferred to another circular orbits with different orbital radius and inclination angle. Figure 6(a) shows the 



violation of the terminal boundary condition according to the number of iterations in a logarithmic scale. 
Red circles denote outer iterations in Newton-Armijo iteration to compute the sensitivity derivatives. The 
error in satisfaction of the terminal boundary condition converges quickly to machine precision after the 
solution is close to the local minimum at around 20th iteration. These convergence results are consistent 
with the quadratic convergence rates expected of Newton methods with accurately computed gradients. 

The neighboring extremal method, also referred to as the shooting method, is numerically efficient in the 
sense that the number of optimiza tion paramet ers is minimized. But, in general, this approach may be prone 
to numerical ill-conditioning (see iBettsI (l20nih 'l. A small change in the initial multiplier can cause highly 
nonlinear behavior of the terminal attitude and angular momentum. It is difficult to compute the gradient for 
Newton iterations accurately, and the numerical error may not converge. However, the numerical examples 
presented in this article show excellent numerical convergence properties. The dynamics of a rigid body 
arises from Hamiltonian mechanics, which have neutral stability, and its adjoint system is also neutrally 
stable. The proposed Lie group variational integrator and the discrete multiplier equations, obtained from 
variations expressed in the Lie algebra, preserve the neutral stability property numerically. Therefore the 
sensitivity derivatives are computed accurately. 



5. Controlled Lagrangian Method for Discrete Lagrangian Systems 



The method of controlled Lagrangians is a procedure for constructing feedback controllers for the stabi- 
lization of relative equilibria. It relies on choosing a parametrized family of controlled Lagrangians whose 
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corresponding Euler-Lagrange flow are equivalent to the closed loop behavior of a Lagrangian system with 
external control forces. The condition that these two systems are equivalent result in matching conditions. 
Since the controlled system can now be viewed as a Lagrangian system with a modified Lagrangian, the 
global stability of the controll ed system can be d etermined directly using Lyapunov stability analysis. 

Th is approach originated inlBloch et al.' ('1997') and was then developed in Auck lv et al. ( 2000( ): Bloch et al 
(|199&.1999 , 2000 , 2001): Hamberg (1999, 2000). A similar approach for Hamihonian controlled systems was 
introduced and further studied in t he work of Blankenste i n, Ortega, va,n der Schaft, Maschke, Spong, and 
their collaborators (see, for example, Maschke et al.l (200ll ): Ortega et al. (2002') and related references). The 
t wo methods we r e shown to be eq uivalen t in IChang et al.l (|2002h and a nonholonomic version was developed 



Zenkov et all (|200d Eioi ). andlBfochl (|20oi l 



Since the method of controlled Lagrangians relies on relating the closed-loop dynamics of a controlled 
system with the Euler-Lagrange flow associated with a m odifled Lagr angian, it is natural to discretize 
this approach through the use of variational integrators. In Bloch et al] |2005l fiooih . a discrete theory of 
controlled Lagrangians was developed for variational integrators, and applied to the feedback stabilization 
of the unstable inverted equilibrium of the pendulum on a cart. 

The pendulum on a cart is an example of an underactuated control problem, which has two degrees of 
freedom, given by the pendulum angle and the cart position. Only the cart position has control forces acting 
on it, and the stabilization of the pendulum has to be achieved indirectly through the coupling between 
the pendulum and the cart. The controlled Lagrangian is obtained by modifying the kinetic energy term, 
a process referred to as kinetic shaping. Similarly, it is possible to modify the potential energy term using 
potential shaping. 

Since the pendulum on a cart model involves both a planar pendulum, and a cart that translates in 
one-dimension, the conflguration space is a cylinder, §^ x M. 

Continuous kinetic shaping: The Lagrangian has the form kinetic minus potential energy 



(50) 



L{qA) 



213(6)9^ 



7s 



(51) 
(52) 



0, 



and the corresponding controlled Euler-Lagrange dynamics is 

d dL dL 

di~d9 ~ 'de 

d dL _ 
dt ds 

where u is the control input. 

Since the potential energy is translationally invariant, i.e., V{q) — V{9), and the relative equilibria 9 = ^g, 
s = const are unstable and given by non-degenerate critical points of V{9). To stabilize the relative equilibria 
9 — 9e, s — const with respect to 9, kinetic shaping is used. The controlled Lagrangian in this case is defined 

by 



(53) 



L^'^iq, q) = Li9, 9, s + T{e)9) + \cj^{T{9)eY 



where t[9) — K.f3{9). This velocity shift corresponds to a new choice of the horizontal space (see |Bloch_et_al 



( 2000t ) for details). The dynamics is just the Euler-Lagrange dynamics for controlled Lagrangian ([5 



(54) 
(55) 



dt 89 



d dL^'' 
dt ds 



= 0, 



= 0. 



Lagrangian ([55)1 satisfies the simplified matching conditions of iBloch et al.l (|200l[ ) when the kinetic energy 
metric coefficient 7 in (j50p is constant. 

Setting u = —d{piT{9)9) /dt defines the control input, makes equations (|52p and ([55)1 identical, and results 
in controlled momentum conservation by dynamics (|5ip and (j52p . Setting a = — I/7K makes equations (|5ip 
and (|54p reduced on the controlled momentum level identical. 

Discrete kinetic shaping: Here, we adopt the following notation: 



9^-1-1/2 



qk + qk+i 



Aqk = qk+i - qk, qk = (Ok, Sk) 
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Then, a second-order accurate discrete Lagrangian is given by, 

Ld{qk,qk+i) = hL{qk+i/2,^qk/h). 
The discrete dynamics is governed by the equations 

dLd{qk,qk+i) dLd{qk-i,qk) _ „ 

dfk + Wk - 

r^ry. dLd{qk,qk+i) , dLd{qk-i,qk) 

(57) + = -life, 

OSk OSk 

where Uk is the control input. Similarly, the discrete controlled Lagrangian is, 

LY{qk,qk+i) = hL^^''{qk+i/2,^qk/h), 

with discrete dynamics given by, 

^^co^ dLY(qk,qk+i) dLY{qk-i,qk) 

^''^ — dfk — + — Wk — 

/ggN dLY{qk,qk+l) ^ dLY{qk-i,qk) ^ ^ 

dsk dsk 

Equation (j59p is equivalent to the discrete controlled momentum conservation: 

Pk = f^, 



where 



Setting 



^ a£;-"(gfc,gfc+i) ^ (l + 7«)/3(gfc+i/2)Agfc+7A5fc 
dsk h 

jA9kT{9k+i/2) - -fAek-iT{9k-i/2) 
Uk = ; 



makes equations ([57]) and ([SS]) identical and allows one to represent the discrete momentum equation (157)) 
as the discrete momentum conservation law pk = p. 

The condition that (I56II57|) are equivalent to ()58II59|) yield the discrete matching conditions. The dynamics 
determined by equations (I56II57P restricted to the momentum level pk — p is equivalent to the dynamics of 
equations ()58ll59p restricted to the momentum level pfc = /i if and only if the matching conditions 

1 p 

a = , ^ = 



JK I + JK 

hold. 

Numerical example: Simulating the behavior of the discrete controlled Lagrangian system involves viewing 
equations (|58ll59p as an implict update map {qk-2,qk-i) '-^ {qk-i,qk)- This presupposes that the initial 
conditions are given in the form [q^, qi); however it is generally preferable to specify the initial conditions as 
(90)9o)- This is achieved by solving the boundary condition, 

dL 

-^(90, go) + DiLdiqo, qi) + -F^ (90, qi) = 0, 

for qi. Once the initial conditions are expressed in the form (gg, qi), the discrete evolution can be obtained 
using the implicit update map. 

We first consider the case of kinetic shaping on a level surface, when k is twice the critical value, and 
without dissipation. Here, h = 0.05 sec, m — 0.14 kg, M = 0.44 kg, and I = 0.215 m. As shown in Figure [71 
the 6 dynamics is stabilized, but since there is no dissipation, the oscillations are sustained. The s dynamics 
exhibits both a drift and oscillations, as potential shaping is necessary to stabilize the translational dynamics. 
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Figure 7. Discrete controlled dynamics with kinetic shaping and without dissipation. The 
discrete controlled system stabilizes the 9 motion about the equilibrium, but the s dynamics 
is not stabilized; since there is no dissipation, the oscillations are sustained. 



6. Future Directions 

Discrete Receding Horizon Optimal Control: The existing work on discrete optimal control has been 
primarily focused on constructing the optimal trajectory in an open loop sense. In practice, model uncertainty 
and actuation errors necessitate the use of feedback control, and it would be interesting to extend the existing 
work on optimal control of discrete systems to the feedback setting by adopting a receding horizon approach. 

Discrete State Estimation: In feedback control, one typically assumes complete knowledge regarding 
the state of the system, an assumption that is often unrealistic in practice. The general problem of state 
estimation in the context of discrete mechanics would rely on good numerical methods for quantifying the 
propagation of uncertainty by solving the Liouville equation. In the setting of H amiltonian syst e ms, t he 
solution of the Liouville equation can be so lved by the method of characteristics ( Scheeres et al.l ( 2007t )V 



This implies that a collocational approach ( Xiu ( 2007[ )) combined with Lie group variational integrators 



and interpolation based on noncommutative harmonic analysis on Lie groups could yield an efficient means 
of propagating uncertainty, and serve as the basis of a discrete state estimation algorithm. 

Forced Symplectic-Energy-Momentum Variational Integrators: One of the motivations for studying the 
control of Lagrangian systems using the method of controlled Lagrangians is that the method provides a 
natural candidate Lyapunov function to study the global stability properties of the controlled system. In 
the discrete theory, this approach is complicated by the fact that the energy of a discrete Lagrangian system 
is not exactly conserved, but rather oscillates in a bounded fashion. 

This can be addressed by considering the symplectic-energy- momentum ( Kane et al. ( 1999l l) analogue of 
the discrete Lagrange-d'Alembert principle, 

Ni N-l 

S'^Ld{qk,qk+i,hk) = ^ [F^ {qk,qk+i,hk) ■ Sqk + F^ [qk,qk+i,hk) ■ 6qk+i], 

k=0 k=0 

where the timestep hk is allowed to vary, and is chosen to satisfy the variational principle. The variations 
in hk yield an Euler-Lagrange equation that reduces to the conservation of discrete energy in the absence 
of external forces. By developing a theory of controlled Lagrangians around a geometric integrator based on 
the symplectic-energy-momentum version of the Lagrange-d'Alembert principle, one would potentially be 
able to use Lyapunov techniques to study the global stability of the resulting numerical control algorithms. 
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